Predicting wildlife corridors for multiple species in an East African ungulate community

Wildlife corridors are typically designed for single species, yet holistic conservation approaches require corridors suitable for multiple species. Modelling habitat linkages for wildlife is based on several modelling steps (each involving multiple choices), and in the case of multi-species corridors, an approach to optimize single species corridors to few or a single functional corridor for multiple species. To model robust corridors for multiple species and simultaneously evaluate the impact of methodological choices, we develop a multi-method approach to delineate corridors that effectively capture movement of multiple wildlife species, while limiting the area required. Using wildlife presence data collected along ground-based line transects between Lake Manyara and Tarangire National Parks, Tanzania, we assessed species-habitat association in both ensemble and stacked species distribution frameworks and used these to estimate linearly and non-linearly scaled landscape resistances for seven ungulate species. We evaluated habitat suitability and least-cost and circuit theory-based connectivity models for each species individually and generated a multi-species corridor. Our results revealed that species-habitat relationships and subsequent corridors differed across species, but the pattern of predicted landscape connectivity across the study area was similar for all seven species regardless of method (circuit theory or least-cost) and scaling of the habitat suitability-based cost surface (linear or non-linear). Stacked species distribution models were highly correlated with the seven species for all model outputs (r = 0.79 to 0.97), while having the greatest overlap with the individual species least-cost corridors (linear model: 61.6%; non-linear model: 60.2%). Zebra was the best single-species proxy for landscape connectivity. Overall, we show that multi-species corridors based on stacked species distribution models achieve relatively low cumulative costs for savanna ungulates as compared to their respective single-species corridors. Given the challenges and costs involved in acquiring data and parameterizing corridor models for multiple species, zebra may act as a suitable proxy species for ungulate corridor conservation in this system.


Introduction
East Africa features an impressive network of protected areas to safeguard landscapes, wildlife assemblages and associated ecosystem processes from human influences [1]. However, many of the protected areas are small, isolated, and not always effective in addressing these goals [2][3][4][5]. As a result, wildlife populations have been declining considerably across the region over the past decades [6][7][8][9]. One of the major concerns regarding the conservation of large mammal populations in Africa and elsewhere is the decline of functional connectivity within ecosystems [10][11][12][13][14][15]. Since connectivity is fundamental for effective landscape-scale conservation of wildlife populations [16], a key strategy to maintain or reverse the loss of functional connectivity is to identify and subsequently protect or restore wildlife corridors-i.e., patches of land that connect two or more protected areas or seasonal ranges of target species [11,17,18].
Identification and delineation of multi-species corridors require spatially-explicit presence, genetic, or movement data for all target species, but such multi-species datasets are rarely available [44][45][46]. Beyond these logistical concerns, modelling multi-species corridors requires the assessment of cross-taxon trade-offs, as it is likely that a corridor designed for one species does not perfectly match corridor requirements for another species [38,47,48]. Failure to account for these cross-taxon differences could compromise the effectiveness of corridors as well as result in additional monetary costs if the primary focus were meant to preserve ideal corridors for all species [38]. Therefore, methods to delineate efficient multi-species corridors that will capture wildlife movement across a mammal community, while limiting the area required to protect them, could provide cost-effective options for connectivity conservation [49][50][51].
In tandem with the understanding of the biological and conservation importance of wildlife corridors, the theory and practice of corridor delineation has made substantial progress over the last two decades [52,53]. Typically, corridor modelling involves (1) the collection of spatially-explicit animal distribution, genetic, or movement data, (2) assembling spatial variables that are hypothesized to be associated with animal presence or movement, (3) fitting appropriate habitat models which allow the prediction of how spatially explicit variables promote or impede animal movement or presence in the landscape, and (4) estimating corridors based on the spatial arrangement of the resistance surface [17,54]. For all these steps, a variety of techniques and data conversion options are available (S1 Fig in S1 File). Because the choices made during each modelling step could have repercussions for the delineation of corridors, comparing the relative effects of different model choices is a crucial step to illustrate uncertainty associated with corridor model parameterization [17,45]. For example, to assess species-habitat associations, a variety of approaches and algorithms are available [55,56]. To address this modelling uncertainty, we harnessed recent advances in species distribution modelling and used ensemble and stacked species distribution models to quantify species-habitat associations [57,58]. Based on the species-habitat associations, scholars estimate resistance surfaces either by the inverse of the habitat suitability (assuming linear relationships between animal movement and landscape resistance) or based on non-linear transformation of the inverse of the squared habitat suitability (assuming non-linear relationships between animal movement and landscape resistance) [59]. As a final methodological choice, the use of either least-cost modelling or circuit theory (the two main approaches in corridor design) may also affect the design of corridors [60,61]. The combined effects of these choices have rarely been quantified [62].
In this paper, we develop a multi-method approach to design corridors for a multi-species assemblage by combing habitat use data from seven ungulate species. Additionally, we assess the value of representative proxy species [63] as a cost-effective method to secure landscape permeability. By following a multi-method approach and quantifying the (relative) consequences of methodological choices for corridor design, we provide an analysis to illustrate uncertainty associated with our corridor modelling approach [45,48]. We hypothesize that (1) species-habitat relationships (and thus movement costs) differ across species, (2) species-specific corridors vary in their spatial configuration, (3) multi-species corridors increase the movement costs for individual species but less so than other species-specific corridors, and (4) choices in the statistical methods and in parameterization of resistance surfaces affect corridor design.

Study area
The study area is located centrally in the Tarangire-Manyara ecosystem, in a 1280 km 2 multiuse area between Lake Manyara (hereafter LMNP) and Tarangire (hereafter TNP) National Parks (centered approximately at 3.5846˚S, 36.0021˚E: Fig 2). The selected study area is ideal for testing multi-species corridor hypotheses because the landscape features a patchwork of protected areas [64] with large mammal populations still occupying land outside fully protected areas [65]. During the wet season, multiple wildlife species leave the national parks (primarily TNP; wildlife in LMNP do not exhibit regular, seasonal movements [2], but do occasionally move in and out of park boundaries [14,66]) and move to areas that are often not formally protected [25, [66][67][68][69][70].
The study area is characterized by a semi-arid climate and the main vegetation type is Acacia-Commiphora savanna [71]. Human population growth and changes in traditional lifestyles have caused substantial expansion of human settlements and, consequently, conversion of natural vegetation to agriculture [72]. The western portion of the study area encompasses parts of the Mto wa Mbu Game Controlled Area, where wildlife falls under the jurisdiction of the Tanzanian Wildlife Authority. Settlements, agriculture, and livestock keeping are technically not permitted [73], yet enforcement is typically weak, and these land uses occur widely in the area [65]. The area also includes Manyara Ranch (hereafter MR), a multiple-use area designated to protect wildlife and support the pastoralist lifestyle of two adjacent villages (Esilalei and Oltukai). MR employs rangers to prevent hunting and ensure that livestock grazing regulations are followed. Livestock herds from adjacent villages are permitted to use the ranch during the dry season whereas ranch livestock graze the area year-round [65]. Settlements, agriculture, and hunting are not permitted on the ranch. The small strip of land between MR and TNP is village lands in the Babati district (MR and Mto wa Mbu Game Controlled Area are situated in the Monduli district) and could not be sampled due to logistical constraints.  (Fig 1). Each transect was surveyed once in April 2015, during the long rainy season. We recorded elephant dung instead of direct sightings as elephants traverse the human-dominated area but are rarely detected directly [74,75]. The study area was divided into a 1 km by 1 km grid; transects bisected each grid north-south, such that parallel transects were separated by 1 km. Along each transect, we recorded sightings of all target species including the GPS coordinates and perpendicular distance from transect using a rangefinder.

Wildlife transects
All research was carried out with permission from the Tanzanian Wildlife Research Institute (TAWIRI) and the Tanzanian Commission for Science and Technology (COSTECH), permit #: 2014-324-ER-2013-191.

Habitat suitability models
To model large mammal species distributions, we created eight raster environmental predictor variables for the study area in ArcMap version 10.7.1 (ESRI, Redlands) (Fig 1). Variables were selected based on hypothesized, biological relationships to the presence of the target species [17]. We calculated the Euclidean distance from bomas (pastoral households and livestock enclosures), LMNP and TNP, and the Manyara River. The locations of bomas were digitized from Google Earth imagery that was captured between 2005 to 2013 [76]. We derived an Enhanced Vegetation Index (EVI) for the study area using the MODIS MOD13Q1 V6 Terra Vegetation Indices 16-Day Global 250 m dataset. Using Google Earth Engine, we took the median EVI value for each 250 m pixel in the study area for all six images captured during the 2015 long rainy season (image dates from April 6 th to May 25 th ). The final four predictor variables were derived from the approximately 1 km resolution FAO Global Land Cover-SHARE database [77]. We used the percentage of density coverage layers for cropland, shrub covered area (shrubland) and tree covered area (woodland). For percentage of density coverage of grassland, we combined the grassland, herbaceous vegetation (aquatic or regularly flooded), and waterbodies layers. We masked out Lake Manyara from the study area using the Global Surface Water Occurrence dataset [78]. As the lake's water level fluctuates over the course of the year, we chose a threshold of 26% water occurrence, which represents the maximum threshold that could be applied without overlapping the transects walked in 2015. We projected all eight raster predictor layers to the Africa Albers Equal Area Conic coordinate system with a 1 km by 1 km cell size. None of the variables had a Pearson correlation coefficient greater than 0.7 (S2 Fig in S1 File), so all were retained for this analysis [79].
To predict species-specific habitat suitability across the study area, we built ensemble species distribution models for each species using the 'SSDM' package [57] in R version 4.0.3 [80]. We sampled environmental data from the eight predictor layers for the GPS locations of direct observations collected along the transects (i.e., presence points). Species with a sufficient sample size (i.e., n > 13) of direct observations [81] included giraffe (n = 34), Grant's gazelle (n = 19), impala (n = 37), Thomson's gazelle (n = 63), wildebeest (n = 36), and zebra (n = 85) (S3 Fig in S1 File). As transects resulted in few direct observations of elephants, we instead used GPS locations of elephant dung piles (n = 803). We performed geographic thinning to account for spatial biases [57], which reduced the number of elephant (n = 696) and Thomson's gazelle (n = 62) observations. Pseudo-absences were generated within the extent of the study area, with number and strategy dependent upon the type of model as recommended in [82]. Ensemble models were built from the highest performing models (AUC � 0.90) of eight different algorithms (generalized linear model, generalized additive model, multivariate adaptive regression splines, generalized boosted regressions model, classification tree analysis, random forest, maximum entropy, artificial neural network, and support vector machines). We evaluated the models using a 70% training/30% evaluation holdout method, with a total of 10 repetitions. To generate a stacked multi-species habitat suitability model, we summed the probabilities of the resulting habitat suitability maps.

Landscape connectivity models
To derive linearly scaled landscape resistance surfaces, we took the inverse of the species habitat suitability predictions from the ensemble models, multiplied them by 100 and added 1 [17,61,83]. This assumes that the cost-weighted distance of travelling across a cell with a predicted habitat suitability of 100% is 1 km (the distance to travel across the cell alone), while the cost-weighted distance of travelling through a cell with a predicted habitat suitability of 1% is equivalent to travelling 100 km through suitable habitat (or 100 times more difficult to cross than a cell with a predicted habitat suitability of 100%). In addition, we squared each of the resulting layers to generate non-linearly scaled cost surfaces [84] to address uncertainty in the relationship between habitat suitability and landscape resistance to movement. For each resulting cost surface layer (linear and non-linear), we modeled landscape connectivity between TNP and LMNP with circuit theory [85] and least-cost methods [86] using Linkage Mapper [87] in ArcMap version 10.7.1 (ESRI, Redlands). As a single cell-wide least-cost path is unlikely to represent wildlife movement, we mapped 10% least-cost corridors using the least-cost output [88,89].

Comparing models
To compare how effective single species are at predicting habitat use and landscape connectivity for the other six species of ungulates in this study, we calculated pairwise Pearson correlations between the species-specific and stacked multi-species habitat suitability, circuit theory and least cost models. We also calculated the percent overlap between each of the single-species predicted least-cost corridor maps.

Habitat suitability models
Ensemble model performance ranged from AUC = 0.91 (elephant) to AUC = 0.94 (Grant's gazelle) (S1 Table in S2 File). Predicted habitat suitability in the study area for all seven species was highest within MR for both the linearly and non-linearly scaled layers (Fig 2 and S4 and S5 Figs in S1 File). The northeastern corner of Lake Manyara was also predicted to be an area of relatively high habitat suitability. The top three variables with the greatest contribution to the models differed for each species (S6 Fig in S1 File and S2 Table in S2 File). Distance to the Manyara River was the predictor variable with the greatest contribution to the habitat suitability models for all species. Distance to bomas had the second highest overall contribution to the SSDM and was one of the top three predictor variables for all species. Distance to a national park had the third highest overall contribution to the SSDM.
Pairwise Pearson correlations between the species-specific and stacked multi-species habitat suitability models indicate a strong positive relationship (>0.7) between the predicted habitat suitability maps across the study area for all species except Thomson's gazelle for both linearly and non-linearly scaled layers (Fig 3 and S4 Table in S2 File). The stacked multi-species model had the highest average pairwise correlation, followed by the zebra model (S3 Table in S2 File).

Landscape connectivity models
The pattern of predicted landscape connectivity across the study area between LMNP and TNP was similar for all seven species regardless of method (circuit theory or least-cost) and scaling of the habitat suitability-based cost surface (linear or non-linear) (Fig 2 and S7-S10 Figs in S1 File). Areas with the lowest landscape resistance ran from the northern tip of TNP across MR, before shifting northwest towards Lake Manyara.
Pairwise Pearson correlations between the species-specific and stacked multi-species landscape connectivity models indicated a strong positive relationship (>0.7) between the predicted circuit theory-derived connectivity maps across the study area for all species except Thomson's gazelle for non-linearly scaled layers (Fig 3 and S3 Table in S2 File). Similarly, the pairwise Pearson correlations showed strong positive relationships between the least-costderived connectivity maps, although the average correlation coefficient was >0.9 for all species and the multi-species model. The stacked multi-species model had the highest average pairwise correlation for both methods, followed by the zebra model (S3 Table in S2 File).

Least-cost corridor models
The landscape connectivity models predicted similar least-cost corridors for each of the seven focal species and the stacked multi-species model (Fig 4 and S11 and S12 Figs in S1 File). Travelling across the study area from Tarangire National Park in the south towards LMNP in the northwest (or vice versa), each species was predicted to move roughly directly across the unprotected lands between TNP and MR. Once in the ranch, each species was predicted to take slightly different optimal movement paths based on the species-specific habitat preferences. Predicted movement patterns appear to maximize the distance travelled through MR. Once exiting the ranch, each of the predicted least-cost corridors roughly converge in the grasslands at the northeastern edge of Lake Manyara and ultimately border the lake before entering LMNP (S11 and S12 Figs in S1 File). Similarly, the optimal multi-species least-cost corridor crosses directly from TNP into MR, maximizes the distance crossed within the ranch, and then directly crosses the Mto wa Mbu Game Controlled Area to reach the northeastern edge of Lake Manyara (Fig 4).
Average pairwise percent overlap between the least-cost corridor maps was generally high (>50%) for all species except elephant, impala (non-linearly scaled only), and Thomson's gazelle (non-linearly scaled only) (Fig 4). The multi-species stacked model had the highest average percent overlap for both the linearly and non-linearly scaled least-cost corridor maps (61.6% and 60.2% respectively), while zebra had the second highest average overlap for the linearly scaled model (59.4%) and giraffe for the non-linearly scaled model (55.0%) (S4 Table in S2 File).

Multi-species corridors
Long-distance movement of vertebrates is a fundamental yet severely threatened process in terrestrial ecosystems across the world [90][91][92]. There is a growing recognition that wildlife corridors should be planned, designed, and implemented for multiple species so that, ideally, the full range of crucial ecosystem processes can be maintained across large spatial scales [36][37][38][39][40][41][42][43]. Conservation planning also ideally allows for movement of a variety of species to facilitate expected range shifts due to anticipated changes in climatic conditions [93][94][95]. Our results suggest that optimizing wildlife corridors by using either a multi-species or a proxy species approach can be a cost-effective method to secure landscape permeability for a large mammal assemblage. Using habitat suitability as a proxy for landscape resistance to movement for seven large mammal species suggests that a single multi-species wildlife corridor would best represent single-species landscape connectivity across the Tarangire-Manyara Ecosystem. However, gathering movement data for multiple species to derive an optimal multi-species corridor might be both time-and cost-prohibitive. Therefore, the identification of a single proxy species that most accurately captures landscape connectivity for an entire or a subset of a species assemblage could reduce these costs, particularly in areas where rapid land conversion is quickly threatening linkages between protected areas.
A previous study identified elephants as a proxy for other large mammal presence and movement in Tanzania [44]. In contrast, we found that zebra might better capture large mammal habitat use and movement in our study area. Differences in the ability of elephants and zebras to predict landscape use of a large mammal community in a typical Miombo ecosystem [44] and an Acacia-Commiphora dominated ecosystem (this study) may be related to actual or field methodrelated differences in mammal community composition in each ecosystem [47]. Irrespective of the underlying mechanism, the poor performance of elephants as a proxy for landscape connectivity in this study emphasizes that conservation planning by proxy species needs a site-specific evaluation to avoid suboptimal conservation outcomes [92]. Landscape connectivity models using both circuit theory and least-cost methods and linearly and non-linearly scaled habitat suitability surfaces for zebra had the highest correlation to connectivity models for the other six species. Furthermore, the percent overlap of the predicted least-cost corridor for zebra compared to the other six species was highest in the linearly scaled model and second highest for the non-linearly scaled model. These findings echo similar results that zebra presence is most closely associated with large mammal species richness in northern Tanzania, while elephants performed poorly as a proxy [47]. The ability of zebra to predict suitable corridors for a range of other herbivore species (including those that differ in feeding strategies) suggest that conservation efforts targeting the protection of corridors for zebra could help conserve landscape connectivity for large mammal assemblages in savanna Africa. Conceptually, space use of habitat generalists may best represent movement of a range of other species. Indeed, space use of zebra usually shows little response to spatial variables, which is indicative of a habitat generalist [96].

Considerations for corridor models
A key assumption in this study is that habitat suitability is analogous to, or at least a reasonable proxy for, landscape resistance. However, species distribution is not necessarily equivalent to movement [54,97]. The relationship between habitat suitability and movement has not been well studied in large mammal species in Africa, but in a recent case study, the distribution of wildebeest was found to be a suitable proxy for their movement [17]. Similarly, the distribution of elephants appears to be highly indicative of their movements as well [23]. However, landscape resistance may differ between day and nighttime [98,99]. As the walking transects were conducted during daylight hours, our study does not capture the potential for differential nighttime habitat use when species might be moving across the landscape to avoid direct encounters with humans and livestock (although this caveat my not apply to elephant space use as the distribution of signs likely captures their actual distribution). Similarly, this study captures large mammal presence during one season of a single year. While the survey was conducted during the rainy season, when dispersal across this landscape generally occurs [66,69,100], spatio-temporal variation in resource availability (e.g., grass, non-alkaline surface water) and the level and extent of the alkaline Lake Manyara might alter landscape resistance. Therefore, optimal wildlife corridors in one year may offer suboptimal paths of movement in following years. However, this study does not primarily aim to offer guidance on where exactly to designate wildlife corridors in the Tarangire-Manyara ecosystem, and so the results from a single season and year are likely valid for testing how best to optimize wildlife corridors for multiple species. However, the consistent importance of distance from nearest boma as a predictor of wildlife presence (S6 Fig) indicates strong spatial avoidance of human settlements by all wildlife species in our study and highlights that human settlements add substantial movement costs to wildlife species [17].
Ideally, landscape resistance and corridor models can be cross-validated with movement data from the seven study species [23,101]. Unfortunately, the only available movement data in the region is for wildebeest [100]. However, none of those five collared individuals entered LMNP, which prevents quantitative comparisons to our modeling efforts. Ideally, and given the importance of landscape connectivity for wildlife species in the fragmented Tarangire-Manyara ecosystem [67,69,70,100,102], additional efforts should be made to capture movement data for a greater number of species to inform corridor conservation in the region.
As we employed a suite of statistical methods and different parameterizations of the resistance surface, we were able to conduct an uncertainty analysis of our modelling choices. While often recommended, such analyses are rarely conducted [48]. Our approach to address uncertainty first considered species distribution model-derived habitat suitability maps as a proxy for landscape resistance. Many methods have been developed to model habitat suitability and rather than pick a single algorithm we chose instead to use an ensemble model in acknowledgment that while all models are flawed, each has their predictive strengths [103]. We then considered uncertainty in the relationship between habitat suitability and landscape resistance to movement. Using both linearly and non-linearly scaled resistance layers, we found that the resulting maps of connectivity were similar when comparing between species with the same species being more (e.g., zebra) or less (e.g., Thomson's gazelle) correlated to others regardless of the shape of transformation. Finally, we considered the impact of connectivity modeling method (least-cost versus circuit theory) on our multi-species comparisons. Similarly, we found that modeling method did not result in major differences in which species better predicted movement for the entire ungulate community.

Conclusions and conservation implications
To model robust corridors for multiple species and evaluate the impact of methodological choices, we developed a multi-method approach and parameterized corridor models for multiple species. Species-habitat relationships and subsequent corridors differed across species, but the pattern of predicted landscape connectivity was similar for all seven species regardless of methodological choices. Stacked species distribution models were correlated with the seven species for all model outputs, while having the greatest overlap with the individual species least-cost corridors.
Connectivity in the fragmented Tarangire-Manyara ecosystem is rapidly declining [17,67,70]. To maintain functional connectivity and current wildlife population sizes (which are already likely below their historic baselines [2,100]) and associated ecosystem services [104] conservation authorities need to implement effective and realistic wildlife corridors. Our fieldand model-based results suggest that either a multi-species corridor or single species (zebra) corridor may effectively facilitate movement of the most abundant herbivore species in this ecosystem while minimizing the land that needs to be set aside for conservation. More generally, this approach highlights the potential feasibility of multi-species corridors for ensuring functional connectivity in savanna ecosystems and emphasizes the need for local evaluation of conservation by proxy approaches.
Supporting information S1 File. Supporting figures. S1 File contains all 12 supporting figures (S1-S12) referenced in the main text along with their associated captions. (DOCX) S2 File. Supporting tables. S2 File contains all four supporting tables (S1-S4) referenced in the main text. Their associated captions follow: S1 Table. Ensemble species distribution model performance metrics. S2 Table. Relative contribution of environmental variables in explaining distribution of wildlife species between Tarangire and Lake Manyara National Parks, northern Tanzania. Some rows may not sum to precisely 100 due to rounding. S3 Table. Pearson correlation coefficients between the predicted habitat suitability maps derived by different methods for modelling the corridor (linear vs. non-linear scaled landscape resistances; least cost path modelling vs. circuit theory). S4 Table. Percent spatial overlap of wildlife corridors modelled with a least-cost algorithm and defining the resistance surface with either a linear or a non-linear scaled landscape resistance value. (XLSX) S1 Data. Supporting data. S1 Data contains the GPS locations of animal sightings used in these analyses.